三维空间点的直线方程拟合(最小二乘法)

最新推荐文章于 2024-05-13 14:52:19 发布

cosmosknight

最新推荐文章于 2024-05-13 14:52:19 发布

阅读量2.5w

点赞数 26

文章标签：经验分享

本文链接：https://blog.csdn.net/qwertyu_1234567/article/details/117918602

版权

本文要解决如下问题：给定一些三维空间点，要求用直线拟合这些点，求出直线方程。

一、直线方程的三种表示方法：

参考：空间直线的点向式方程 - 百度文库 (baidu.com)

空间直线方程 - 百度文库 (baidu.com)

1.一般式：

它实际上表示，直线是两个平面的交线，因此可以由两个平面方程得到，即：

$\left\{\begin{matrix} A_{1} X + B_{1} Y + C_{1} Z + D_{1} = 0\\ A_{2} X + B_{2} Y + C_{2} Z + D_{2} = 0 \end{matrix}\right.$

2.点向式（标准方程）：

$\frac{x-x_{0}}{m} = \frac{y-y_{0}}{n}=\frac{z-z_{0}}{p}$

(m, n, p) 为直线方程的方向向量；(x0, y0, z0) 为直线上的一个点。需要注意的是(x-x0, y-y0, z-z0)的方向和方向向量是平行的，也因此推导出了上面的方程。

3.参数方程：

$\frac{x-x_{0}}{m} = \frac{y-y_{0}}{n}=\frac{z-z_{0}}{p} = t$

由此就可以得到：

$\left\{\begin{matrix} x = x_{0} + mt \\y=y_{0} + nt \\ z = z_{0} +pt\end{matrix}\right.$

二、利用最小二乘，由三维空间数据点拟合直线方程：

参考文献：《三维空间点中基于最小二乘法的分段直线拟合方法》薛丽红，2015年7月，齐齐哈尔学报，第31卷第4期

相应链接：https://www.doc88.com/p-8189740853644.html

联立上面的方程，就可以得到：

公式中，划红线的2,是该篇文章的笔误，应该是n, n表示数据点的个数。

由此我们就可以得到k1, b1, k2, b2。如果想要把它转换成空间直线标准方程，可以先任意取一个z0值，根据k1,b1,k2,b2 就可以求出x0, y0, 然后任意取一个p值，根据k1, k2, 就可以求出m, n值，代入到标准方程即可(所以，标准方程不是唯一的，需要从标准方程的定义去理解为什么不是唯一的)。

三、python 代码实现：

##  由空间3维点拟合出一条直线
def linear_fitting_3D_points(points):
    '''
    用直线拟合三维空间数据点。
    
    参考https://www.doc88.com/p-8189740853644.html  
    《三维空间点中基于最小二乘法的分段直线拟合方法》 薛丽红，2015年7月，齐齐哈尔学报，第31卷第4期
    注意; 文中的公式推导有误，k1,b1,k2,b2中的系数2， 应该为n，n表示数据点的个数。
    
    直线方程可以转化成如下形式（具体见上面的文献）：
    x = k1 * z + b1
    y = k2 * z + b2
    
    Input:
        points    ---   List， 三维空间数据点，例如：
                        [[2,3,48],[4,5,50],[5,7,51]]
                    
    返回值是公式系数 k1, b1, k2, b2
    '''

    #表示矩阵中的值
    Sum_X=0.0
    Sum_Y=0.0
    Sum_Z=0.0
    Sum_XZ=0.0
    Sum_YZ=0.0
    Sum_Z2=0.0

    for i in range(0,len(points)):
        xi=points[i][0]
        yi=points[i][1]
        zi=points[i][2]

        Sum_X = Sum_X + xi
        Sum_Y = Sum_Y + yi
        Sum_Z = Sum_Z + zi
        Sum_XZ = Sum_XZ + xi*zi
        Sum_YZ = Sum_YZ + yi*zi
        Sum_Z2 = Sum_Z2 + zi**2

    n = len(points) # 点数
    den = n*Sum_Z2 - Sum_Z * Sum_Z # 公式分母
    k1 = (n*Sum_XZ - Sum_X * Sum_Z)/ den
    b1 = (Sum_X - k1 * Sum_Z)/n
    k2 = (n*Sum_YZ - Sum_Y * Sum_Z)/ den
    b2 = (Sum_Y - k2 * Sum_Z)/n
    
    return k1, b1, k2, b2

如有错误，请指正。谢谢！